library(cyclestreets)
library(mapview)
library(pct)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(stats19)
## Data provided under OGL v3.0. Cite the source and link to:
## www.nationalarchives.gov.uk/doc/open-government-licence/version/3/
library(stplanr)
library(tidyverse)
## -- Attaching packages ------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.2.1 v purrr 0.3.2
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0
## -- Conflicts ---------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(tmap)
library(tmaptools)
library(dplyr)
library(leaflet)
library(classInt)
library(ggplot2)
region_name = "liverpool-city-region"
max_distance = 7
lines_test = get_pct_lines(region_name)
zones_all = get_pct_zones(region_name)
L_original_lsoa = get_pct_lines("liverpool-city-region", geography = "lsoa")
Liverpool_filter = filter(L_original_lsoa, lad_name1 %in% c("Liverpool")) # filter out Liverpool only
tmap_mode("view")
## tmap mode set to interactive viewing
#3a Bicycle Dependent Desire Lines
bicycle = Liverpool_filter %>%
mutate('Percentage Cycling' = (bicycle) / all * 100) %>%
filter(e_dist_km < max_distance)
tm_shape(bicycle) +
tm_lines("Percentage Cycling", palette = "RdYlBu", lwd = "bicycle", scale = 9)
## Legend for line widths not available in view mode.
#3b Car Dependent Desire Lines
car_dependent = Liverpool_filter %>%
mutate(`Percent Drive` = (car_driver) / all * 100) %>%
filter(e_dist_km < max_distance)
tm_shape(car_dependent) +
tm_lines("Percent Drive", palette = "-RdYlBu", lwd = "all", scale = 9)+
tm_layout(legend.text.size = 1.2)
## Legend for line widths not available in view mode.
#note: these desire lines can be created for any travel mode, adapt the code above to show desired mode
Top10_Liverpool = Liverpool_filter %>% top_n(10, bicycle)
#filter the columns
Top10_Liverpool_data = Top10_Liverpool %>% select(id, geo_code1, geo_code2, geo_name1, geo_name2, lad11cd1, lad11cd2, lad_name1, lad_name2, all, bicycle, foot, car_driver, car_passenger, motorbike, train_tube, bus, taxi_other, govtarget_slc, dutch_slc)
#you will require a cyclestreets API to run this code
#code = Top10_Liverpool_route = stplanr::line2route(Top10_Liverpool_data, route_fun = stplanr::route_cyclestreet)
Top10_Liverpool_route <- st_read("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/Data/Top10_liverpool_route.shp")
## Reading layer `Top10_liverpool_route' from data source `C:\Users\Holly.Mizser-Jones\OneDrive - Arup\UCL\CASA005\Assessment\Data\Top10_liverpool_route.shp' using driver `ESRI Shapefile'
## Simple feature collection with 11 features and 20 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: -2.99466 ymin: 53.34071 xmax: -2.83573 ymax: 53.41017
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
#join the Government Target, Go Dutch and Existing Bicycle trips to the raw route data
Top10_Liverpool_route$govtarget = Top10_Liverpool_data$govtarget_slc
Top10_Liverpool_route$dutch = Top10_Liverpool_data$dutch_slc
Top10_Liverpool_route$bicycle <- Top10_Liverpool_data$bicycle
#map graph showing the variation in hilliness / distance between the Government Target and Go Dutch scenarios
#difference between gradient and distance
distances = 1:20
hilliness = 0:5
uptake_df = data.frame(
distances = rep(distances, 6),
hilliness = rep(hilliness, each = 20)
)
p_govtarget = uptake_pct_govtarget(
distance = uptake_df$distances,
gradient = uptake_df$hilliness
)
p_godutch = uptake_pct_godutch(
distance = uptake_df$distances,
gradient = uptake_df$hilliness
)
uptake_df = rbind(
cbind(uptake_df, scenario = "govtarget", pcycle = p_govtarget),
cbind(uptake_df, scenario = "godutch", pcycle = p_godutch)
)
library(ggplot2)
ggplot(uptake_df) +
geom_line(aes(
distances,
pcycle,
linetype = scenario,
colour = as.character(hilliness)
)) +
scale_color_discrete("Gradient (%)")+
labs(x = "Distance (km)", y = "Propensity to Cycle (Multiplier)")
#read in LSOA shapefile for UK with IMD 2019 data joined
Liverpool_lsoa_shp <- st_read("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/Shp/Indices_of_Multiple_Deprivation_IMD_2019/Indices_of_Multiple_Deprivation_IMD_2019.shp")
## Reading layer `Indices_of_Multiple_Deprivation_IMD_2019' from data source `C:\Users\Holly.Mizser-Jones\OneDrive - Arup\UCL\CASA005\Assessment\Shp\Indices_of_Multiple_Deprivation_IMD_2019\Indices_of_Multiple_Deprivation_IMD_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32844 features and 66 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -6.418524 ymin: 49.86474 xmax: 1.762942 ymax: 55.81107
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
#filter for Liverpool
Liverpool_filter_shp = filter(Liverpool_lsoa_shp, LADnm %in% c("Liverpool"))
#visualise the routes on top of the IMD polygon data
library(leaflet)
pal <- colorFactor(get_brewer_pal("-YlGnBu", n = 10), domain = Liverpool_filter_shp$IMD_Decile, n = 10)
pal2 <- colorFactor(get_brewer_pal("-YlGnBu", n = 10), domain = Liverpool_filter_shp$HDDDec, n = 10)
IMDandroutes<- leaflet() %>%
# add basemap options
addProviderTiles(providers$Esri.WorldGrayCanvas, group = "Esri Grey") %>%
addTiles(group = "OSM") %>%
#add our Borough polygons, linking to the tables we just made
addPolygons(data=Liverpool_filter_shp,
color="white",
weight = 2,
opacity = 1,
dashArray = "3",
fillOpacity = 0.6,
fillColor = ~pal(Liverpool_filter_shp$IMD_Decile),
group = "Index of Multiple Deprivation")%>%
# add a legend for boroughs
addLegend(pal = pal,
values = Liverpool_filter_shp$IMD_Decile,
group=c("Index of Multiple Deprivation"),
title ="IMD Decile",
position ="bottomleft")%>%
#add our ward polygons, linking to the tables we just made
addPolygons(data=Liverpool_filter_shp,
color="white",
weight = 2,
opacity = 1,
dashArray = "3",
fillOpacity = 0.6,
fillColor = ~pal2(Liverpool_filter_shp$HDDDec),
group = "Health Deprivation and Disability Decile")%>%
# add a legend for wards
addLegend(pal = pal2,
values = Liverpool_filter_shp$HDDDec,
group=c("Health Deprivation and Disability Decile"),
position ="bottomleft",
title ="HDD Decile")%>%
#add routes
addPolylines(data=Top10_Liverpool_route, color = "black",
weight = 2.5, opacity = 1, group = c("Top 10 Cycle Routes in Liverpool")) %>%
# specify layers control
addLayersControl(
baseGroups = c("Esri Grey Canvas", "OSM"),
overlayGroups = c("Index of Multiple Deprivation", "Health Deprivation and Disability Decile","Top 10 Cycle Routes in Liverpool"),
options = layersControlOptions(collapsed = FALSE))%>%
hideGroup(c("Index of Multiple Deprivation", "Health, Deprivation and Disability Decile", "Top 10 Cycle Routes in Liverpool"))
# show us the map
IMDandroutes
7a Health Intersect
#intersect the Top Routes with the IMD polygon data
IMD_route_intersect <- st_join(Top10_Liverpool_route, Liverpool_filter_shp, join = st_intersects)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
#export this to csv and use Microsoft Excel to calculate the average for each route
write.csv(IMD_route_intersect,"C:/Users/Holly.Mizser-Jones/Documents/UCL/CASA005/Assessment//IMD_route_average.csv", row.names = FALSE)
#read the IMD average back into R
IMD_route_average <- read.csv("C:/Users/Holly.Mizser-Jones/Documents/UCL/CASA005/Assessment/IMD_average.csv")
#merge with the route data
route_merge <-merge (Top10_Liverpool_route,
IMD_route_average,
by.x="length",
by.y="length",
no.dups = TRUE)
7b Calculate the change in Marginal METS (Goodman et al., 2019)
#the change in mMETS was calculated in excel, the below code joins the resultant increases to the routes
#link MET data to Top_10_routes
mMET <- read.csv("C:/Users/Holly.Mizser-Jones/OneDrive - Arup/UCL/CASA005/Assessment/liverpool_mMET.csv")
#merge mMET with Routes
route_mMET <-merge (route_merge,
mMET,
by.x="Route",
by.y="Route",
no.dups = TRUE)